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Simple and Efficient Numerical Evaluation 
of Near-Hypersingular Integrals 


Patrick W. Fink, Member, IEEE, Donald R. Wilton, Fellow, IEEE, Michael A. Khayat, Member, IEEE 


Abstract — Simple and efficient numerical procedures for 
evaluating the gradient of Newton-type potentials are presented. 
Convergences of both normal and tangential components of the 
gradient are examined. The convergence of the vector potential 
is also examined, and it is shown that the scheme for handling 
near-hypersingular integrals also is effective for the nearly 
singular potential terms. 

Index Terms — Numerical integration, Integral equations, 
Computation theory . 

i. Introduction 

R ecently, significant progress has been made in the 
handling of singular and nearly- singular potential 
integrals that commonly arise in the Boundary Element 
Method (BEM) [1-3]. To facilitate object-oriented 
programming and handling of higher order basis functions, 
cancellation techniques are favored over techniques involving 
singularity subtraction. However, gradients of the Newton- 
type potentials, which produce hypersingular kernels, are also 
frequently required in BEM formulations. As is the case with 
the potentials, treatment of the near-hypersingular integrals 
has proven more challenging than treatment of the limiting 
case in which the observation point approaches the surface. 
Historically, numerical evaluation of these near- 
hypersingularities has often involved a two-step procedure: a 
singularity subtraction to reduce the order of the singularity, 
followed by a boundary contour integral. Departure from 
object-oriented philosophy and heightened complexity in the 
handling of higher order basis functions are even more 
problematic in this context. Thus, there is ample need for 
cancellation-type techniques for efficient numerical evaluation 
of the gradient of the potential. 

Progress in the development of efficient cancellation-type 
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procedures for the gradient potentials was recently presented 
[2, 4]. The approach in deriving these procedures closely 
parallels that which led to the so-called “Arcsinh” 
transformation [1]. That is, to the extent possible, a change of 
variables is chosen such that the Jacobian of the 
transformation cancels the singularity. However, since the 
gradient kernel involves singularities of different orders, we 
extend the approach to require that the transformation reduces 
the order of the singularity and leaves remaining terms that are 
analytic. The terms “normal” and “tangential” are used herein 
with reference to the source element. 

For the singular potential terms, the integral is defined and 
well-behaved when the observation point is contained within 
the source element. However, for the gradient potential terms 
we must be content with numerical evaluation in which the 
observation point remains some finite, but arbitrarily small, 
distance, d, above the source element. Procedures for 
extracting the normal component of the gradient in the limit as 
the observation point approaches the source element are well 
known and are not discussed herein. 

The work in [2] presents the desired singularity order 
reduction for the gradient terms and enables extremely 
efficient computation of the normal gradient component as the 
distance z tends toward zero. However, convergence as 
presented in [2] is strongly dependent upon z. Also, the 
tangential gradient components are not addressed in [2], 
primarily because the techniques therein do not permit 
noticeably efficient evaluation. Also, since computational 
formulations often involve the numerical evaluation of the 
potentials in addition to their gradients, it is highly desirable 
that a single integration procedure efficiently handle both. 
Techniques to extend and improve the methods in [2] were 
recently presented [4] and are the primary subject of this 
paper. 

II. Integrals with near-Hypersingularities 
A. Theory 

In this section, the transformations previously presented [2, 
4] are reviewed. Further modifications to improve 
convergence, of both the normal and tangential gradient 
components, are established. Cartesian and spherical 
coordinates are used herein to simplify notation and to readily 
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convey mathematical behavior. In practice, use of simplex 
coordinates defined on the parent element is often preferred, 
and conversion of quadrature points on sub-elements follows 
simple procedures previously outlined [1], 

In order to simplify the presentation and evaluation of the 
proposed methods, we consider the following integral 
representative of the actual integrals with respect to the 
convergence of numerical integration: 
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where R = |r - r'| is the distance between an observation point 

at r and source points on triangular domains V . Here, the 
function b is intended to capture the behavior of the basis 
function. Following previous approaches [1-5], a nearby 
observation point is projected onto the surface containing the 
source parent triangular element. The parent triangle is then 
split into three subtriangles about the projected observation 
point. The geometry of a typical subtriangle and its local 
rectangular coordinate system with origin at the projected 
observation point is shown in Fig. 1. 

The general form of the transformed integral over the 
subtriangle is given by 

p* prc / x 
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where g{x, y) is a product of the basis and the Green’s 


function. Ideally, the transform’s Jacobian J would exactly 
cancel the singular (static) part of the kernel for constant 
source densities, and the subtriangle would map into a 
rectangular domain such that repeated Gauss-Legendre 
integration of low order may be used. A summary of two 
transformations previously investigated [2, 4] is given in 
Table 1, followed by a synopsis of their properties. 



Note that both transformations result in terms, p , which arise 
in the linear basis function components and in the x- and y- 
components of the vector, R. These terms, which are non- 
analytic in v, must be handled carefully, and there are other 
detractions associated with the R ’ transformation to be noted. 
It is extremely efficient for evaluating the static normal 
component of the gradient (k = 0). However, it is very poor 
in capturing any of the other terms [2]; that is, p and higher 
powers of R , due to the near-singularities introduced in these 
terms by the transformation. Moreover, we present here a 
means, using the R 2 transformation, to handle the normal 
component of the gradient with the same high efficiency, but 


none of the detractions, characterizing the R 3 method. For this 
reason, the R 3 is not considered further in this paper. The R 2 
transformation results in smooth functions, but can result in 
slow convergence when the ratio of the end boundary value of 
v to the initial boundary value of v is large. 

TABLE 1: 

SUMMARY OF TRANSFORMATIONS 



TRANSFORMATION 

J(u,v) 

INTEGRATION LIMITS 

Radial-Angular 

u = (f> 

R 2 

6, n 

u, . = In tan — — 

L .t 2 

R 2 
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Radial-Angular 

U = (f> 

-R i 

U LU = Au 

R 3 

\z\ 
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TT 
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Jz + (A/sin u) 


Radial-Angular R 2 : 

• Introduces R 2 in the Jacobian to reduce singularities in R 

• R becomes a smooth function of the integration variable, 
v 

• Introduces terms, p , that are non-analytic at R = z 
Radial- Angular R ’: 

• Introduces R 3 in the Jacobian to reduce singularities in R 

• Exactly integrates the normal gradient component of the 
static potential with one sample point per subtriangle 

• Introduces nearly-singular components, (|z|/v) m . 

• Introduces terms, p , that are non-analytic at R = z 


To improve the convergence using the R 2 transformation, 
we consider the introduction of a disk of radius p centered at 

the projection point as shown in Figure 2, along with the three 
subtriangles. Since the transformations listed in Table 1 
already employ radial and angular boundaries, inclusion of the 
disk element is seamless. The introduction of this fourth sub- 
element provides three important advantages. First, the 
azimuthal variation of the gradient integrand can be integrated 
exactly with only three sample points, assuming linear basis 
functions (completeness order p = 0), or with two sample 
points for constant basis functions. The 3-point azimuthal rule 
comprises three equally-spaced angles and equal weights - a 
Riemann sum; e.g., {(/ } = 2zr {l/ 6, 1 / 2 , 5 / 6 } and 

{w } = { 1 / 3 } . It exactly integrates the set of functions 
{f.} = sin (j>, cos <!>, sin 2<j>, cos 2 <£} over the circle. For 
potentials, only two sample points are required for exact 
integration of the azimuthal component of the integrand with 
linear bases, or one sample point for constant bases. The 
second key advantage pertains to the non-analyticity of p . 
Since the azimuthal components are integrated exactly, all odd 
powers of p , which are non-analytic and are multiplied by 
odd powers of the sine or cosine functions, are entirely 
eliminated over the disk. 
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Figure 2. Disk sub-element centered at projection point. 

The third key advantage pertains to the radial variation of 
the integrand, and we now show that accurate integration of 
the radial terms over the disk can also be achieved with great 
efficiency. Upon application of the Radial-Angular R 2 
transformation, the integrand of (1) can be written as 

‘g = ) ( l + J kR ) e jkR ^ > (3) 

in which dependencies upon the variables of integration have 
been suppressed. Since all terms with odd powers of the 
cylindrical radius, p , are eliminated upon integration over 
the disk, the integrand is reduced to factors with the following 
powers of R: 

& = , 1, R,} {l - jkR - k 2 R 2 - ...} , (4) 

where the complex exponential has been expanded. We now 
consider integration of a single term of power m and make use 
of a normalized variable of integration, A , 

V u 1 

I m =\e dv = A\e m ^ +Vl) dA, (5) 

v, 0 

where 

m = - 1,0,..; v tu = In (z,R 0 ),R 0 = ^ p\ + z 2 , 
and A = v u - v ( . 

Equation (5) can then be rewritten 

i 

/, = Ajz”'e mAIna dA (6) 

0 

where 

« - Rjz . (7) 

We now consider means to numerically integrate the dominant 
terms in (4) efficiently. We assume that the disk radius is 
small relative to a wavelength, and thus address only the 
constant and linear components in the expansion of the phase 
term. The set of ^-dependent terms is then reduced to 

f -in (a)A . in (a)A 2\n(a)A . ) 

= \e y , 1, e ' , e v M. (8) 

For efficient integration of these terms over the disk, we 
employ a 2-point quadrature rule, which integrates the terms 
in (8) exactly, providing that (7) is satisfied. Since (7) 
pertains to a specific ratio of radius to z, we create a table of 


2-point quadrature rules to cover the desired range of as. 
Each table entry corresponds to a specific value of a , and we 
choose the entry, a , nearest to the desired value given by (7). 
We then modify the disk radius to match the selected table 
entry, 

R 0 = z a, (9) 

so that all terms represented in (8) are integrated exactly over 
the disk. 

If the table has a sufficient number of entries, the modified 
disk radius is sufficiently close to the desired radius that there 
are no adverse consequences, as is demonstrated in the next 

section. For a table accommodating values 1 < a < 10 5 , a 

maximum ratio ( z ) = 2 , and a deviation from the 

desired disk radius below ten percent, the total table size is 
103 quadrature rules. 

It remains to determine the desired value of disk radius. 
Since integration over the disk is particularly efficient, and 
integration external to the disk is prone to effects of the non- 
analyticity of p at R = z, larger values of disk radius appear 
advantageous. Extension of the disk region beyond the 
subtriangle boundaries (as shown in Fig. 2) is readily 
accommodated and requires no additional considerations. In 
this case, outlying regions of the disk acquire a negative 
weight due to the reversal of integration between the circular 
boundary and the outer subtriangle edge. Initial trials, 
however, reveal no significant advantages to extending the 
radius beyond the triangle boundaries, and it is convenient to 
set the disk radius 

p =h , (10) 

r 0 min v 7 

the minimum subtriangle height (Figure 1). 

B. Results 

The triangle shown in Figure 3 is used as a test case for 
three different heights, z, of the observation point: 10' 6 m, 0.01 
m, and 0.2 m. For all but one case, as noted below, the disk 
radius equals 0.1 m. The wavelength is 10.0 m so that the 
longest edge is 0. 1 wavelengths. The basis function variation 
is chosen as b = x + y (Figure 3), and the R 2 transformation is 
used for all trials. 

Except as noted below, the specialized quadrature schemes 
presented above for the u- and v-directions are applied for 
integration over the disk region. The disk radius is set to 0. 1 
m, except for one case to demonstrate insensitivity to this 
parameter. Over the truncated triangular regions, Gauss- 
Legendre is used to integrate the transformed integrals along 
both radial and angular directions. Each truncated subtriangle 
was analyzed individually to determine near-optimal sampling 
ratios along the two directions (m,v) in order to achieve a near 
monotonic convergence for both the normal and tangential 
gradient components. Convergence trends of the normal and 
tangential components are shown in Figures 4 and 5, 
respectively, with number of sample points indicating the 
combined total for disk and all three subtriangle regions. 

As seen in Figure 4, accurate integration of the normal 
gradient component is extremely efficient as z tends to zero. 
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The crossed square symbols in this figure denote disk 
quadrature schemes with fewer than six sample points. In the 
extreme case, just 1 sample point (in the disk) is used for the 
whole source triangle. In contrast to previous work 
demonstrating great accuracy with 1-2 points [2], the 
introduction of the disk permits continued convergence of the 
normal component with more sample points, as well as 
convergence of the tangential component. The cross symbols 
in Figure 4 represent use of a disk radius p 0 = 0.09, and it 

can be seen that the effects are minimal. As z becomes 
sufficiently large, the v-components can be integrated 
efficiently over the disk with a two-point Gauss-Legendre 
scheme. For example, in Figures 4 and 5, for the case z = 0.2 
and the range of sample points shown, very little difference 
was observed when the v-components were integrated with a 
two-point Gauss-Legendre scheme compared to the two-point 
specialized quadrature described in the preceding section. 



Figure 3. Triangle geometry for trials. 


Figure 6 shows the convergence results for the potential 
using the same quadrature rules applied to the gradient. The 
crossed square symbols indicate trials in which the 3x2 
quadrature rule over the disk is replaced by a 3x3 rule in order 
to assess eventual limitations due to the truncation of the 
phase term implicit in (8). It should be noted that one could 
improve the convergence for the potential by employing the 
same concepts (i.e., specialized quadratures over the disk) and 
optimizing specifically for the potential. For example, only 
two points are required for the w-integration over the disk, and 
the term with the negative exponent in (8) would be replaced 
by the next positive integer exponent. Further work is 
required to compare this approach with another recent 
advance for integration of 1/R-type singularities [3]. 
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Figure 5. Convergence of tangential gradient component. 



Figure 6. Convergence of the potential using the gradient- 
based quadrature rules. 
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Figure 4. Convergence of normal gradient component. 
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